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3PI ALGORITHM FOR SPIRAL CT 

This invention claims the benefit of priority to U.S. Provisional Application, 
Serial No. 60/430,802 filed, December 4, 2002, and is a Continuation-In-Part of United 
States Patent Application Serial No. 10/389,534 filed March 14, 2003 which is a 
Continuation-In-Part of Serial No. 10/389,090 filed March Continuation-In-Part of Serial 
No. 10/143,160 filed May 10, 2002 now U.S. Patent 6,574,299, which claims the benefit 
of priority to U.S. Provisional Application 60/312,827 filed August 16, 2001. 

FIELD OF INVENTION 
This invention relates to computer tomography, and in particular to processes, 

methods and systems for reconstructing three-dimensional images from the data obtained 

by spiral scans using 3PI algorithm. 

BACKGROUND AND PRIOR ART 
Over the last thirty years, computer tomography (CT) has gone from image 

reconstruction based on scanning in a slice-by-slice process to spiral scanning. From the 

1970s to 1980s the slice-by-slice scanning was used. In this mode the incremental 

motions of the patient on the table through the gantry and the gantry rotations were 

performed one after another. Since the patient was stationary during the gantry rotations, 

the trajectory of the x-ray source around the patient was circular. Pre-selected slices 

through the patient have been reconstructed using the data obtained by such circular 

scans. From the mid 1980s to present day, spiral type scanning has become the preferred 

process for data collection in CT. Under spiral scanning a table with the patient 

continuously moves through the gantry that is continuously rotating about the table. At 

first, spiral scanning has used one-dimensional detectors, which receive data in one 
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dimension (a single row of detectors). Later, two-dimensional detectors, where multiple 
rows (two or more rows) of detectors sit next to one another, have been introduced. In 
CT there have been significant problems for image reconstruction especially for two- 
dimensional detectors. In what follows the data provided by the two-dimensional 
5 detectors will be referred to as cone-beam (CB) data or CB projections. 

For three-dimensional (also known as volumetric) image reconstruction from the 
data provided by a spiral scan with two-dimensional detectors, there are two known 
groups of algorithms: Exact algorithms and Approximate algorithms, that each have 
known problems. Under ideal circumstances, exact algorithms can provide a replication 

10 of an exact image. Thus, one should expect that exact algorithms would produce images 
of good quality even under non-ideal (that is, realistic) circumstances. However, exact 
algorithms can be known to take many hours to provide an image reconstruction, and can 
take up great amounts of computer power when being used. These algorithms can require 
keeping considerable amounts of cone beam projections in memory. Additionally, some 

1 5 exact algorithms can require large detector arrays to be operable and can have limits on 
the size of the patient being scanned. 

Approximate algorithms possess a filtered back projection (FBP) structure, so 
they can produce an image very efficiently and using less computing power than Exact 
algorithms. However, even under the ideal circumstances they produce an approximate 

20 image that may be similar to but still different from the exact image. In particular, 

Approximate algorithms can create artifacts, which are false features in an image. Under 
certain circumstances these artifacts could be quite severe. 

The first group of algorithms that are theoretically exact and have the shift- 
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invariant FBP structure was disclosed in U. S. Patent Application 10/143,160 filed May 
10 3 2002, now U.S. Patent 6,574,299, which claims the benefit of priority to U.S. 
Provisional Application 60/312,827 filed August 16, 2001. A shortcoming of these 
algorithms is that they operate in the 1PI mode. This means that if the detector array is 
5 large in the axial direction, then one has to translate the patient through the gantry very 
quickly in order to use all of the detector. However, rapid motion of the patient is very 
impractical for obvious reasons. On the other hand, if the patient moves slowly, only part 
of the detector is used. This leads to undesirable consequences: part of the x-ray dose is 
wasted, discretization artifacts are enhanced, noise stability is reduced, etc. 

10 SUMMARY OF THE INVENTION 

A primary objective of the invention is to provide 3PI algorithms for 
reconstructing images of objects that have been scanned in a spiral fashion with two- 
dimensional detectors. For image reconstruction at any given voxel these algorithms 
require a longer section of the spiral than the 1PI algorithms of U. S. Patent Application 

15 10/143,160 filed May 10, 2002, now U.S. Patent 6,574,299, which is incorporated by 

reference, which claims the benefit of priority to U.S. Provisional Application 60/312,827 
filed August 16, 2001 . Consequently, the new algorithms allow to slow the patient down 
by about a factor of three, but still use the same size detector array. 

A first preferred embodiment of the invention uses a five overall step process for 

20 reconstructing the image of an object under a spiral scan. In a first step a current CB 

projection is measured. Next, families of lines are identified on a detector according to a 
novel algorithm. Next, a computation of derivatives between neighboring projections 
occurs and is followed by a convolution of the derivatives with a filter along lines from 
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the selected families of line. Next, using the filtered data, the image is updated by 
performing back projection. Finally, the preceding steps are repeated for each CB 
projection until an entire object has been scanned. This embodiment works with keeping 
several (approximately 2-4) CB projections in memory at a time and uses one family of 
5 lines. 

For the second embodiment, different families of lines can be used in combination 
with keeping several CB projections in memory. 

Modifications of these embodiments are possible, that will allow keeping only one 
CB projection in computer memory at a time. This can be done analogously to what was 
10 done in U. S. Patent Application 10/143,160 filed May 10, 2002, which is incorporated by 
reference, which claims the benefit of priority to U.S. Provisional Application 60/3 12,827 
filed August 16, 2001. 

Consequently, the new algorithms allow to slow the patient down by about a 
factor of three, but still use the same size detector array. 
1 5 Further objects and advantages of this invention will be apparent from the 

following detailed description of the presently preferred embodiments, which are 
illustrated schematically in the accompanying drawings. 

BRIEF DESCRIPTION OF THE FIGURES 
20 Fig. 1 shows a typical arrangement of a patient on a table that moves within a rotating 
gantry having an x-ray tube source and a detector array, where cone beam projection data 
sets are received by the x-ray detector, and an image reconstruction process takes place in 
a computer with a display for the reconstructed image. 
Fig. 2 shows an overview of the basic process steps of the invention. 
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Fig. 3 shows stereographic projection of the spiral onto the detector plane 
Fig. 4 shows the detector plane with various projections and important lines 
Fig. 5 shows the boundary curve 

Fig. 6 shows the continuous illumination case, x/R = (0,0.25,0) 
5 Fig. 7 shows the interrupted illumination case, x/R = (-0.5,0,0) 
Fig. 8 shows the points where various critical events occur 

Fig. 9 shows the distribution of weights inside the 5IP domain in the case of continuous 
illumination 

Fig. 10 shows the distribution of weights inside the 5IP domains in the case of interrupted 
10 illumination 

Fig. 1 1 shows the family of filtering lines parallel to L 0 

Fig. 12 shows the family of filtering lines tangent to T ±1 

Fig. 13 shows how to choose filtering lines depending on the location of x 

Fig. 14 shows the family of filtering lines tangent to T ±2 

15 Fig. 15 shows the filtering lines and the associated constants c m in different cases: x e F^ 
(top left panel), xeF 2 (top right panel), x e F 3 and above T, (middle panel), xeG x 
(bottom left panel), xeG 2 (bottom right panel). 
Fig. 16 shows possible locations of points s }9 s 29 s 3 

Fig. 17 shows the filtering lines and the associated constants c m in different cases: xsF x 

20 (left panel), x <e F 2 (right panel). 

Fig. 18 is a three substep flow chart for finding families of lines for filtering, which 
corresponds to step 20 of Fig. 2. 
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Fig. 19 is a seven substep flow chart for preparing for filtering, which corresponds to step 
30 of Fig. 2. 

Fig. 20 is a seven substep flow chart for filtering, which corresponds to step 40 of Fig. 2. 
Fig. 21 is a five substep flow chart for the first part of back-projection, which corresponds 
5 to step 50 of Fig. 2. 

Fig. 22 is a three substep flow chart for the second part of back-projection, which 
corresponds to step 50 of Fig. 2. 

DESCRIPTION OF THE PREFERRED EMBODIMENTS 
Before explaining the disclosed embodiments of the present invention in detail it 

10 is to be understood that the invention is not limited in its application to the details of the 

particular arrangements shown since the invention is capable of other embodiments. 

Also, the terminology used herein is for the purpose of description and not of limitation. 

This invention is a continuation-in-part of U. S. Patent Application 10/143,160 

filed May 10, 2002, now U.S. Patent 6,574,299, which is incorporated by reference, 

15 which claims the benefit of priority to U.S. Provisional Application 60/3 12,827 filed 

August 16, 2001. 

The invention will now be described in more detail. 

THEORETICAL BACKGROUND 
First we introduce the necessary notations. Let 

20 C={ye R 3 :y } = Rcos(s\y 2 = Rsin(s) 9 y 3 = s{hI2n\s e R}, (1) 

where h > 0 be a spiral, and U be an open set strictly inside the spiral: 

Ucz{xeR 3 :x;+xl<r 2 }, 0<r<R, (2) 

5 2 is the unit sphere in R 3 , and 
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J3(s,x):= 



\x~y(s)\ 



,xeU,seR, n(x,<f) := {y e R 3 :(y-x) = 0}, (4) 



that is D f (y, fi) is the CB transform of / . Given (x, £)eUx(R 3 , 0) , let 
^ = s j (£>€'X),j = l 9 2 9 ..., denote points of intersection of the plane n(x, £) with C . 
5 Also, y(s) := . For /? e S 2 , /? x denotes the great circle {a e S 2 :a /3 = 0} . 

Fix any x e R 3 , where / needs to be computed. In order to compute / (x) we will use a 
section of the spiral of finite extent, which is to be determined later. For now it will be 
denoted C(x) . The corresponding parametric interval is denoted /(x) . The main 
assumption about C(x) is the following. 



Property CI. (Completeness condition) Any plane through x intersects C(x) at 
least at one point. 

An important ingredient in the construction of the inversion formula is weight 
function n(s,x,a\a e fi L (s,x) . The function n can be understood as follows, x and a 
15 determine the plane n(x,a) , and the weight n assigned to y(s) e Yl(x,a) n C depends 
on the location of x . In view of this interpretation we assume n(s,x,a) = n(s,x,-a) . The 
main assumptions about n are the following. 
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Property Wl. Normalization condition: 



20 



^ n(Sj , x, a) = 1 a.e.on S 2 ; 



(5) 



j 
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Property W2. There exist finitely many C 1 functions 
a k (s,x) € /3 L (s,x),se I(x) , such that n(s,x,a) is locally constant in a neighborhood of 
any (s, a), where s e I(x) and a e p l (s,x),a g Kj k a k (s,x) . 
Denote 

<f>(s,x,6):=sgn(ay{s))n(s,x,a), a = a(s,0)e 0 L (s,x), (6) 
Under assumptions CI, Wl, and W2 the following general inversion formula is derived 
in the paper by A. Katsevich "A general scheme for constructing inversion algorithms for 
cone beam CT," International Journal of Mathematics and Mathematical Sciences, Vol. 
21, pp. 1305-1321 (2003) : 

/(x) = — Lf y^i^L 

x t ^-D f (y(q),cosr/3(s,x) + sinra L (s,x,0 m ))\ -^-ds, (7) 
* oq '^sin^ 

where 6 m e[0,n) are the points where <j>(s,x,0) is discontinuous, and c m {s,x) are values 
of the jumps: 

c m (s,x) := \im(t(s t x,e m +£)-t(s,x,e m -e)). (8) 

£:->0 

Unless n is chosen appropriately, the inversion formula is not necessarily of the FBP 
type. 

3PI LINES AND THEIR PROPERTIES 
Suppose that the x-ray source is fixed at y(s 0 ) for some s 0 g R . Since the 

detector array rotates together with the source, the detector plane depends on s 0 and is 

denoted DP(s Q ). It is assumed that DP(s 0 ) is parallel to the axis of the spiral and is 

tangent to the cylinder y] +y\ = R 2 (cf. (1)) at the point opposite to the source. Thus, the 
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distance between y(s 0 ) and the detector plane is 2R (see Fig. 3). Introduce coordinates 

in the detector plane as follows. Let the d x -axis be perpendicular to the axis of the spiral, 

the d 2 -axis be parallel to it, and the origin coincide with the projection of y(s Q ) . Project 

stereographically the upper and lower turns of the spiral onto the detector plane as shown 
5 in Fig. 3. This gives the following parametric curves: 
, sin(s-.s 0 ) ... h s-s 0 

d x (s) = 2R- v , d 2 (s) = — -± (9) 

1-cos(j-5 0 ) n 1-cos(j-j 0 ) 

p + 2n(j-l)<s-s 0 <2nj^p,j>l, or (10) 

10 p + 2nj<s-s 0 <2;r(y + l)-p,y<-l, (11) 

where p is determined by the radius of support of the patient: p = 2cos~ l (r/R) (cf. (2)). 

These curves are denoted Y p j = ±1,±2,... (see Fig. 4). x denotes the projection of x . 

Connect an arbitrary source position y(s Q ) to all points y(s) on the spiral, where 

27r<s-s Q <47r or -An < s-s 0 < -In. (12) 

15 This gives two surfaces, which are denoted Sy PI (s 0 ) and S 3 L PI (s 0 ) . The region bounded 
by the two surfaces and the cylinder x, 2 +x 2 2 = R 2 is denoted V 3PI (s 0 ) . Let x be fixed. If 
s 0 is sufficiently small, then Slf 1 (s 0 ) is below x . As s 0 increases, two cases are 
possible. In the first one, known as continuous illumination, x enters V 3PI (s Q ) through 
S 3PI (s 0 ) and leaves V 3PI (s Q ) through Sl PI (s 0 ) . Clearly, the above procedure yields a 

20 unique 3PI interval [b Q (x\t 0 (x)] and the corresponding unique 3PI line L 3PI (x) . In the 
second case, known as interrupted illumination, x enters and leaves V 3P1 (s Q ) several 
times. More precisely, x intersects each of the surfaces Sl PI (s 0 ) and S 3 L PJ (s 0 ) exactly 
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three times. Therefore, the above procedure now gives three 3PI lines L 3PI (x\i = 1,2,3 . 
The corresponding values of the parameter are denoted &,(jc),f,(x),z = 1,2,3 . Due to the 
symmetry of the spiral, if x leaves (enters) V 3PI (s 0 ) when s 0 =Z>,(x), then x enters 
(leaves) V* PI (s 0 ) when s 0 = r. (x), / = 1,2,3. Consider the plane jc 3 = 0 . The boundary 
between the two cases is shown in Fig. 5. We see that the curve has no self intersections, 
so it divides the open disk x 2 x + x\ < 1 into two regions: X x and X 3 . In the central one, 
denoted X x , there is one 3PI line for each x . In the exterior one, denoted X 3 , there are 
three 3 PI line for each x . 

The following properties can be established. If there is only one 3PI line for a 
point x , then 

x e V 3P1 (s Q ) « s 0 e I™ (x) := [Vol (13a) 
2;r</,.-6 ; <4;r, i = 1,2,3. (14a) 
If there are three 3 PI lines for a point x , then 

xGK 3F/ ^ 0 )«5 0 G/ 3W (x):=[6 p 6Ju[&3,rJu[r 2 ,r 3 ^ (13b) 
2n<t i -b i <47r, i = 1,2,3. (14b) 
In terms of the detector plane equations (13a) and (13b) imply that x is between 
T 2 and T_ 2 if and only if s e I 3PI (x) . Therefore, using the analogy with the 1PI case, we 
can call this region on the detector the 3 PI window. In the 1PI case, the parametric 
interval bounded by the endpoints of the 1PI line of x is called the 1PI parametric interval 
of* . Similarly, I 3PI (x) is called the 3PI parametric interval of x . 
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AN AUXILIARY CONSTRUCTION 
Suppose first that continuous illumination takes place. Consider the following two 

curves on the surface of the unit sphere S 2 . This first curve consists of all unit vectors a 

orthogonal to L 0 3PI (x) and is denoted by A . The other curve consists of vectors 

5 a(s)=± ±z2^m. Je/ .» (x)i (l5) 

\(x-y(s))xy(s)\ 

and is denoted by T . It is not convenient to represent these curves directly on S 2 , so they 
will be shown on the plane using spherical coordinates (0 P 0 2 ) defined by 

S 2 3<z(s) = (cos <9, sin ^ 2J sin ^ sin 6> 2 , cos ^ 2 ), -n<G x <7r,O<0 2 <n. (16) 
10 Since both a and -a define the same plane, we can restrict 0 ] to any interval of length 
n and "glue" the opposite boundaries using the identification 

W^S^+^U^-^). (17) 
A typical situation for x/R = (0,0.25,0) is shown in Fig. 6. It is very convenient to think 

15 about points on S 2 not only as unit vectors, but also as planes. Each aeS 2 corresponds 
to a unique plane through x with normal vector a . This correspondence is one-to-one if 
vectors with opposite orientation are identified. 

Together A and T split S 2 into several domains: D„Z) 2 ,.... By construction, for a fixed 
i , the number of points in C 3PI (x) n Tl(x, a) is the same for all a € D i , . If £> contains it 
20 intersection points (IPs), it will be called a £ IP domain. 

Suppose now that interrupted illumination takes place. In a similar fashion, define 
several curves on the surface of S 2 . The first three curves are obtained by considering 
unit vectors perpendicular to each of the three 3PI axes. They are denoted A k ,k = 1,2,3 . 

The second set of curves is obtained by restricting s in (15) to the intervals [b x ,b 2 ] , 

11 
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[b 3 ,t x ] , and [t 2 ,t 3 ] • These curves are denoted T n , T M , and T 23 , respectively. A typical 
situation for x/R = (-0.5, 0,0) is shown in Fig. 7. 

It is clear that if x, and x 2 are close to each other, they will share similar 
diagrams. By this we mean that by a smooth sequence of transformations one diagram can 
be converted into the other, and in corresponding domains the number of IPs and their 
distribution over the subintervals in I 3PI stays the same. An essential change is possible 
only in a neighborhood of x where a "critical event" occurs: three boundaries intersect 
each other at one point on 5 2 . These points can be found numerically and the results are 
shown in Fig. 8. The smallest distance from any of these points to the center of rotation is 
r0« 0.618# . Thus, in situations of interest in medical applications (r FOV < 0.5R ), only 

the following two cases can happen: continuous illumination as shown in Fig. 6 and 
interrupted illumination as shown in Fig. 7. 

CONSTRUCTION OF THE WEIGHT FUNCTION n 

For any s e P PI (x) determine the weight function n(s,x,a\a e f3 x (s,x) , according to 
the following rule: 

• In IIP domains the only IP gets weight n = 1 ; 

• In 3 IP domains each IP gets weight n = 1/3 ; 

• In 5IP domains two IPs get weight n = 0 , and all the remaining IPs get weight 
« = l/3 (see Figs. 9 and 10). 
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FIRST RECONSTRUCTION ALGORITHM 
Once the weights are found, we have all the ingredients needed for constructing 

the algorithm. First of all, for image reconstruction at x we take the 3PI spiral segment of 
x : C(x) = C 3PI (x) and I(x) = P P1 (x) (cf. Section 1). As follows from (6), (7), for each 
s € P PI (x) the filtering directions are determined by finding the discontinuities of 
0(s 9 x y a) = sgn(a • y(s))n(s,x,a) . The study of these discontinuities leads us to the 
following three families of lines. 

The first family consists of lines parallel to the spiral tangent. This family is 
denoted £ 0 . The lines and the associated coefficients c m are shown in Fig. 1 1 . Obviously, 
for any x there is only one line Le£ 0 that contains x . 

In Fig. 12 we see the second family of filtering lines. It consists of lines tangent to 
T ±1 and is denoted £ x . Here for each x within the 1PI window we take two lines from 
£ x . These lines are determined from the rule explained in Fig. 13. This figure also shows 
the associated coefficients c m . 

In Fig. 14 (left panel) we see the family of filtering lines tangent to T ±2 . This 
family is be denoted £ 2 . Now, for any x e F } u F 2 u F 4 u F 5 there might be more than 
one line L e L 2 that contains x . The unique line is determined from the following rule 
(see Fig. 14, right panel). If x e F x u F 4 , the point of tangency is to the right of x . If 
xeF 2 KjF 5 , the point of tangency is to the left of x . In all cases c m = 1/3 . 

Fig. 15 summarizes the information contained in Figs. 1 1-14. It shows all possible 
cases where x might be and all the associated filtering directions and constants c m . In all 
cases the direction of filtering as assumed to be from left to right. 
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It follows from our construction (cf. Fig. 15) that 6 m (s,x) and c m (s,x) in the general 
reconstruction formula (7), (8) depend on x only via fi(s,x) . Therefore, we can replace 
x by J3(s,x) in the arguments of c m and a x and rewrite (7) in the form 

m = ^l»^X^^\ g c m (s^ m (s^x))ds^ (18) 

V m (s,0):= rl-^/CX^cos^ + sin^ 1 ^^^))! -$L, (i 9 ) 

x, 2 +x 2 2 < (0.618if) 2 . (20) 
Step 3 1. Fix a line L from the said set of lines obtained in Step 20. 
Note also that given a filtering line L , all x whose projections belong to L and satisfy 

the rules mentioned above share the same filtering line L (cf. Fig. 15). Thus, (19) 

becomes a convolution, (18) becomes backprojection, and the algorithm (18), (19) is of 

the convolution-based FBP type. This can be seen similarly to U. S. Patent Application 

10/143,160 filed May 10, 2002, which is incorporated by reference. Other methods and 

techniques for backprojection can be used. See additionally, for example, US Patent 

6,574,299 to Katsevich, which is incorporated by reference. 

IMPROVED RECONSTRUCTION ALGORITHM 
Let y/(t) be any smooth function defined on R and with the properties ^(0) = 0, 

y/\t) > 0,t e R . Define the new family of lines J? 2 ' by requesting that any given line 

L g J? 2 ' has three points of intersection with T ±I u T ±2 : s l9 s 29 s^ , and these points 
satisfy: 

S } -S=y/(S 3 -S 2 \ S + 27T <S 3 <S + 47T, (21) 

s 3 -s 2 =if/(s x -s), s-An < s z < s-2n. (22) 
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Recall that s is the current source position. The lines L e £ 2 can be parameterized, for 
example, by s^ln <| s 3 1< An . Location of s x and s 2 depends on where s 3 is and is 
illustrated in Fig. 16. 

Using the properties of y/ it is easy to establish that for each x e F x kj F 2 u F 4 u F 5 

there is a unique L e £ 2 that contains jc . Also, if x L 2 , then s 2 ,s 3 -> 2 A and s x s . 
Similarly, if x L c _ r 2 , then s 2 ,s 3 -2 A and s { 5 . 

More importantly, in the case x e F x u F 2 u F 4 u F 5 it is possible to reduce the 
number of filtering directions from two to one. The additional benefit is the improved 
detector usage. Thus, the top two panels shown in Fig. 15 should be replaced by the 
diagrams shown in Fig. 17. The case when x appears below L 0 (i.e., jc e F 4 u F 5 ) can be 
obtained from Fig. 17 by symmetry with respect to the origin in the detector plane. 

The form of the inversion formula (18), (19), remains the same. The first 
difference is that M(s, (5) = 1 (instead of 2) when jc 6 F x kj F 2 kj F 4 u F 5 . The second 
difference is the direction of filtering, which is now determined from (21), (22). 

GENERAL OVERVIEW OF THE RECONSTRUCTION ALGORITHMS 
Fig. 2 shows an overview of the basic process steps 10, 20, 30, 40, 50 of the 

invention. The steps will now be described. 

Step 10, Load the current CB (cone beam) projection into computer memory. Suppose 
that the mid point of the CB projections currently stored in memory is> > ( l y 0 ) . 
Step 20. Finding families of lines for filtering. 
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Using either Fig. 15 or, additionally Fig. 17 in case of the improved algorithm, identify 
the required families of lines. Then, choose a discrete set of lines from each family. 
Step 30. Preparing for filtering. 

Parameterize points on the said lines selected in Step 20 by polar angle. Using 
interpolation compute the derivative of the CB data (dl dq)D f {y(q\P) \ q=SQ at points /? 
on the said lines that correspond to discrete values of the polar angle. 
Step 40. Filtering. For each line identified in Step 20 convolve the data for that line 
computed in Step 30 with filter 1 /s'my . 

Step 50. Back-projection. For each reconstruction point x back-project the filtered 
data found in Step 40 according to equation (18). Then go to Step 10, unless there are no 
new CB projections to process or image reconstruction at all the required points x have 
been completed. 

Now we describe the algorithm in detail following the five steps 10-50 shown in Fig. 2. 

Step 10. Load the current CB (cone beam) projection into computer memory. Suppose 
that the mid point of the CB projections currently stored in memory isy(s 0 ) . The detector 
plane corresponding to the x-ray source located at y(s 0 ) is denoted DP(s 0 ) . 

Step 20. Finding families of lines for filtering. 

The set of lines can be selected by the following substeps 21, 22, and 23. 

Step 21. From the family of lines J? 0 choose an equidistant set of lines that are 
parallel to the spiral tangent and that cover the projection of the region of interest 
onto the detector plane located between T 2 and T_ 2 (see Fig. 11). 
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Step 22. From the family of lines £ x choose a discrete set of lines that are tangent 
to T, and T_j (see Fig. 12). The extreme left point of tangency on r, does not 
have to extend beyond the projection of the region of interest onto the detector 
plane. Similarly, the extreme right point of tangency on r_, does not have to 
extend beyond the projection of the region of interest onto the detector plane. 
Step 23. From the family of lines £ 2 choose a discrete set of lines that are tangent 
to T 2 and r„ 2 (see Fig. 14, left panel). In both cases the points of tangency do 

not have to extend beyond the projection of the region of interest onto the detector 
plane. 

In case of the improved algorithm, instead of the lines tangent to T 2 and T_ 2 , we 
choose a discrete (say, equidistant) set of values for s 3 on the curves T 2 and 

T_ 2 and then determine the lines L e J? 2 ' by solving equations (21), (22). On both 
curves the points s 3 do not have to extend beyond the projection of the region of 
interest onto the detector plane. 

Step 30, Preparing for filtering 

Step 3 1 . Fix a line L from the said set of lines obtained in Step 20. 

Step 32. Parameterize points on the said line by polar angle y in the plane 

through y(s 0 ) and L . 

Step 33. Choose a discrete set of equidistant values y. that will be used later for 
discrete filtering in Step 40. 
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Step 34. For each y } find the unit vector (3j which points from y(s Q ) towards the 
point on L that corresponds to y. . 

Step 35. Using the CB projection data D f (y(q),&) for a few values of q close to 

s 0 find numerically the derivative (d I dq)D f (y{q),&) \ q=SQ for all 0 = fi. . 

Step 36. Store the computed values of the derivative in computer memory. 
Step 37. Repeat Steps 31-36 for all lines L identified in Step 20. This way we 
will create the processed CB data corresponding to the x-ray source located at 

Step 40, Filtering 

Step 41. Fix a line L from one of the families of lines £ m identified in Step 20. 

Step 42. Compute FFT of the values of the said processed CB data computed in 

Step 30 along the said line. 

Step 43. Compute FFT of the filter \/s'my 

Step 44. Multiply FFT of the filter 1/sin/ (the result of Steps 43) and FFT of the 

values of the said processed CB data (the result of Steps 42). 

Step 45. Take the inverse FFT of the result of Step 44. 

Step 46. Store the result of Step 45 in computer memory. 

Step 47. Repeat Steps 41-46 for all lines in the said families of lines. This will 

give the filtered CB data *¥ m (s Q9 /3j) , where m stands for the line family number 

from which L was selected, m = 0,1, 2 . 
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By itself the filtering step is well known in the field and can be implemented, for 
example, as shown and described in U.S. Patent 5,881,123 to Tarn, which is 
incorporated by reference. 

Step 50. Back-projection 

Step 51. Fix a reconstruction point x , which represents a point inside the patient 
where it is required to reconstruct the image. 

Step 52. If s 0 belongs to I 3PI (x) , then the said filtered CB data affects the image 
at x and one performs Steps 53-58. If s 0 is not inside the interval P PI (jc) , then the 
said filtered CB data is not used for image reconstruction at jc. In this case go 
back to Step 5 1 and choose another reconstruction point. 
Step 53. Find the projection x of x onto the detector plane DP(s Q ) and the unit 
vector /?O 0 ,x), which points from y(s 0 ) towards x. 

Step 54. Using Figs. 1 1, 13, and the right panel of Fig. 14 identify the lines from 
the said families of lines and points on the said lines that are close to the said 
projection x . If x is above 1% or below Zf r 2 , then in the case of the improved 
algorithm use equations (21), (22) and Fig. 16 to find the filtering lines close to 
the said projection x . This will give a few values of *¥ m (s O9 0j) for /?. close to 

Step 55. With interpolation estimate the value of m (s 0 ,j3(s 0 ,x)) fr° m the said 
values of *F m ($ 0 ,/? y ) for fi. close to 0(s o ,x). 
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Step 56. Compute the contribution from the said filtered CB data to the image 
being reconstructed at the point x by multiplying *F w (s 0 ,/?(j 0 ,x)) by 

-c m (s,/3(s 09 x))/[47r 2 1 x-y(s 0 ) |] . The appropriate backprojection coefficient 

c m is selected using Figs. 1 1, 13, and the right panel of Fig. 14 (see also Fig. 15 for 

the summary). If x is above L c 2 r or below LZ 2 , then in the case of the improved 

algorithm use Fig. 17 to find the appropriate backprojection coefficient c m . 

Step 57. Add the said contributions to the image being reconstructed at the point 

x according to a pre-selected scheme (for example, the Trapezoidal scheme) for 

approximate evaluation of the integral in equation (18). 

Step 58. Go to Step 51 and choose a different reconstruction point x . If all 

reconstruction points x have been processed, go to Step 59. 

Step 59. Go to Step 10 and load the next CB projection into computer memory. 

The image can be displayed at all reconstruction points x for which the 
image reconstruction process has been completed (that is, all the subsequent CB 
projections are not needed for reconstructing the image at those points). Discard 
from the computer memory all the CB projections that are not needed for image 
reconstruction at points where the image reconstruction process has not 
completed. The algorithm concludes when the scan is finished or the image 
reconstruction process has completed at all the required points. 

For example, if the detector does not provide all the data which is required 
for the 3 PI algorithm, one can employ various techniques for estimating the 
missing data. In this case by the detector (respectively, cone beam projection) we 
mean the virtual detector (respectively, virtual cone beam projection), which 
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includes both measured and estimated data. If one is able to estimate the missing 
data exactly, then the 3PI algorithm will produce exact reconstruction. In this 
sense we still talk about exact reconstruction under real circumstances, when 
missing data are found approximately. 

5 

While the invention has been described, disclosed, illustrated and shown in 
various terms of certain embodiments or modifications which it has presumed in practice, 
the scope of the invention is not intended to be, nor should it be deemed to be, limited 
thereby and such other modifications or embodiments as may be suggested by the 
10 teachings herein are particularly reserved especially as they fall within the breadth and 
scope of the claims here appended. 
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